Fluctuation statistics of mesoscopic Bose-Einstein condensate: reconciling the master 
equation with the partition function to revisit the Uhlenbeck-Einstein dilemma 
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The atom fluctuations statistics of an ideal, mesoscopic, Bose-Einstein condensate is investigated 
from several different perspectives. By generalizing the grand canonical analysis (applied to the 
canonical ensemble problem), we obtain a self-consistent equation for the mean condensate parti- 
cle number that coincides with the microscopic result calculated from the laser master equation 
approach. For the case of a harmonic trap, we obtain an analytic expression for the condensate 
particle number that is very accurate at all temperatures, when compared with numerical canon- 
ical ensemble results. Applying a similar generalized grand canonical treatment to the variance, 
we obtain an accurate result only below the critical temperature. Analytic results are found for 
all higher moments of the fluctuation distribution by employing the stochastic path integral for- 
malism, with excellent accuracy. We further discuss a hybrid treatment, which weds the master 
equation/stochastic path integral analysis with the results obtained based on canonical ensemble 
quasiparticle formalism [V. V. Kocharovsky et al., Phys. Rev. A 61, 053606 (2000)], producing 
essentially perfect agreement with numerical simulation at all temperatures. 

PACS numbers: 03.75.Hh,05.30. Jp 



I. INTRODUCTION 

The problem of fluctuation statistics in mesoscopic 
Bose-Einstein condensate (BEC) is a rich one - even, 
somewhat surprisingly, for the noninteracting Bose gas, 
as well as for the interacting Bose gas . Take, for exam- 
ple, a trapped gas of bosons, and lower the temper- 
ature below the critical temperature Tc to form a con- 
densate. Count the number of atoms that are in the 
condensate. Do this many times at the same tempera- 
ture and total boson number N to get statistics. What is 
the distribution of the condensed atoms, and how does it 
change when moving the temperature through Tel The 
noninteracting, mesoscopic, Bose gas is a problem whose 
fluctuation characteristics are still not completely under- 
stood. Contrary to standard lore, we stress that the fluc- 
tuations are not Gaussian, even in the thermodynamic 
limit. The discussion goes clear back to the Uhlenbeck- 
Einstein dilemma 0, Q : Uhlenbeck's criticism that Ein- 
stein's expression for the average boson number (no) had 
an abrupt cusp at Tc- The manner in which this cusp 
is smoothed by fluctuations is of great interest for meso- 
scopic systems, which is one focus of the present paper. 
Even more subtle is the question of the higher order mo- 
ments of the condensate fluctuations, especially in the 
vicinity of Tc- 

Despite the conceptual simplicity of the above ques- 
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tion, the fluctuation statistics are not known analytically, 
because while the canonical ensemble partition function 
can be formally written, it can be accurately calculated 
only numerically. We note that the work by Holthaus and 
Kalinowski obtains accurate approximations for the first 
few moments, employing a refined saddle point approxi- 
mation of the canonical partition function However, 
the equation for the saddle point still has to be solved 
numerically. 

The time-dependent master equation is an excellent 
alternative approach to attack this question, developed 
in the pioneering papers of Scully b\ (hereafter referred 
to as CNBl) and Kocharovsky, Scully, Zhu and Zubairy 
(hereafter referred to as CNB2). This approach is 
capable of investigating dynamical phenomena, as well 
as being a viable option when we include interactions, 
or when the system is out of equilibrium. The master 
equation approach, within the canonical framework, has 
successfully calculated moments of the fiuctuation dis- 
tribution to good accuracy. Other advantages of this 
method are that it provides a simple physical picture of 
the condensation process with a gas of cold atoms via a 
detailed balance of the elementary processes of heating 
and cooling, and that the theory shows a vivid analogy 
between atom BEC and the photon laser by demonstrat- 
ing that the BEC and laser master equations are identi- 
cal. This analogy can explain the atom-laser linewidth 
Q, and might also be able to explain the transient dy- 
namics of the phase transition to BEC. In the master 
equation approach, the cooling of the BEC is accom- 
plished by the non-condensate interaction with a thermal 
environment. In typical BEC experiments, the cooling is 
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FIG. 1: (Color online) (Left). The average boson number in the condensate is plotted versus temperature. The Master equation 
approach (solid line) is compared with numerical simulation of canonical results (dots), as well as the thermodynamic limit 
(dashed line). (Right). The variance of the distribution is plotted versus temperature. The grand canonical result (dashed 
line) deviates near and below the critical temperature, leading to the "grand canonical catastrophe". The result of ter Haar 
[To| (dash dot) and Politzer flU (small dots) (derived in the thermodynamic limit) are also shown. The results are shown for 
A'^ = 200 particles in a harmonic trap. 



accomplished with evaporative techniques. However, in 
the experiments by Reppy and co-workers the BEC 
was cooled by a thermal reservoir associated with the Vy- 
cor glass host - exactly the situation we presently con- 
sider. The master equation takes a given bath model 
for the calculation, so one natural question is if a differ- 
ent bath model alters the physics of the BEC statistics. 
While different models will certainly alter the dynamics 
of the condensation, we expect the statistics to be robust, 
just as the analogous partition function is. 

How does all of this relate to traditional statistical 
physics, in particular the microcanonical, canonical, and 
grand canonical formalism? Recall that Einstein pre- 
dicted BEC by analyzing the grand canonical formula 
. Both grand canonical ensemble and canonical ensem- 
ble approaches describe the occurrence of BEC equally 
well in the dramatic rise of the mean condensate num- 
ber (no) in the ground state. However, near the critical 
temperature and below, the grand canonical approach 
gives grossly different results for the variance and higher 
moments, the so-called "grand canonical catastrophe" il- 
lustrated in Fig. njj. 

The question of fluctuations in Bose-Einstein conden- 
sation has attracted renewed interest. Grossmann and 
Holthaus 0, ^3 previously obtained an analytical ex- 
pression for the mean condensate particle number in 
a harmonic trap for finite N using the grand canoni- 
cal approach. However, their analytical expression does 



not show a smooth crossover near Tc- Grossmann and 
Holthaus and Ketterle and Druten found the 
crossover behavior only by numerical solution of the 
equation for the condensate particle number. Politzer 
also obtains an approximate expression for the variance 
in the thermodynamic limit 11]. 

The purpose of the present paper is to report advances 
on several fronts. One goal is to connect the master 
equation formalism with the partition function analysis 
in thermal equilibrium. To this end, we review the mas- 
ter equation formalism, and derive an approximate result 
for the mean condensate particle number, that can be ex- 
pressed in terms of elementary functions. This result is 
then connected with standard statistical mechanics by 
extending the grand canonical approach (applied to the 
canonical ensemble problem!) to obtain an approximate 
result for the mean condensate particle number that is 
quite accurate through the critical temperature. This is 
done by deriving a self-consistent equation for the con- 
densate particle number by eliminating the chemical po- 
tential. The equation may be approximately solved, and 
the solution exactly coincides with the approximate mas- 
ter equation result. 

Next, we investigate fluctuations, a more subtle prob- 
lem. The analogous treatment for the variance via the 
generalized grand canonical analysis works quite well 
for low temperatures, but fails for temperatures near or 
greater than the critical temperature. In order to obtain 
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a more accurate analytic result from the master equation 
approach, we first reformulated the master equation as 
a stochastic_^ath integral. The stochastic path integral 
formalism 15, allows for an accurate approximation 
of all cumulants of the fluctuation statistics in terms of 
elementary functions, by applying the saddle-point ap- 
proximation, whose large parameter is the number of 
condensed atoms. It is straightforward to generalize the 
results, and present the fluctuation statistics of a gen- 
eral boson master equation, with arbitrary heating and 
cooling coefficients. 

Finally, by combining the master equation/stochastic 
path integral predictions with the previous results of V. 
V. Kocharovsky et al. based on canonical ensemble quasi- 
particle formalism ^Y?\ (hereafter referred to as CNB3), 
which is very accurate at low temperature, we are able 
to formulate a hybrid theory, whose results are in near- 
perfect agreement with exact numerical simulation. The 
idea is to use the low-T results of CNB3 to fix the values 
of three master equation parameters via the first three 
moments. Then the master equation predictions for sev- 
eral higher moments and cumulants are compared with 
the numerical results at all temperatures. They are in 
near-perfect agreement. 

The paper is organized as follows. In Sec. II we review 
the master equation approach, and derive an approxi- 
mate solution for the mean atom number in the conden- 
sate. This solution is connected with standard partition 
function methods in Sec. Ill by extending the grand 
canonical analysis, and deriving a self-consistent equa- 
tion for the mean atom number, which coincides with 
the result of Sec. II. In Sec. IV, we consider a harmonic 
trap, and give improved analytic results for the heating 
coefficient. We go on to apply the extended grand canon- 
ical analysis to the variance of the distribution. In Sec. 
V, we give the stochastic path integral treatment of the 
problem, which yields approximate expression for all cu- 
mulants in terms of the master equation parameters. In 
Sec. VI, we introduce a hybrid theory, combining the pre- 
dictions of the master equation/stochastic path integral, 
with those of CNB3, to obtain excellent agreement with 
numerical simulation. We conclude in Sec. VII. Appen- 
dices A and B contain the details of the calculations for 
the grand canonical approach, and appendix C contains 
the details of the stochastic path integral calculations. 



cesses of a Bose-Einstein condensate is given by |5| 

+ HnoUoPno - Hno + l{no + l)pno + l], (1) 

where in the low temperature limit 

Kno^N-no, Hno='H, (2) 



fe>0 

and K is a rate constant. The steady state solution for 
the distribution is 

~ (N~no)\e'^r{N+l,n)' ^ ' 

where 

r{N + l,x) = / t'^e-^dt (5) 

J X 

is an incomplete gamma-function. Higher order moments 
may be found via (ng) = J2no ^oPno- We introduce the 
notation /is = ((no — "■o)") for the central moments of 
the distribution. The equation for the probability distri- 
bution can also be converted to an equation for the mean 
atom number (no) to find 

(no) = ~n[{-N + 1 + 7i)(no) - iV + (n^)]. (6) 

Approximating (np) ~ {no)'^, and considering the sta- 
tionary case, the solution for (no) is 

(no) ~ i (iv - w - 1 + ^{N -n-i)^ + m^ . (7) 

In contrast to CNBl, this approximate answer can be 
simply expressed in terms of H and iV, yet still catches 
the smooth transition in the vicinity of Tc (where N — 
Ti. ~ VN), see Fig. Elbelow. We will show in the next sec- 
tion that this same (approximate) result can be derived 
within the context of the grand canonical formalism, and 
later in Sec. V, how to derive similar results for all higher 
moments from the stochastic path integral formalism. 



II. MASTER EQUATION APPROACH 

Let us review the results of the master equation ap- 
proach. Consider a trapped Bose gas, whose single- 
particle energy levels are e^. The master equation for- 
malism begins with the elementary transitions of atoms 
between the condensate and non-condensate. We are pri- 
marily interested in the condensate statistics, and there- 
fore focus on the probability Pn^ of having no atoms in 
the condensate, given that there are N total atoms. The 
master equation describing the heating and cooling pro- 



Ill. GRAND CANONICAL TREATMENT 

In order to see how the same result for the average con- 
densate particle number emerges from a modified grand 
canonical treatment, we recall that the grand canonical 
formula for the average of the total particle number N is 

fc=0 fc=0 
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where the sum runs over all states of energy e^, l3~^ = 
ksT and the chemical potential is related to the mean 
condensate particle number (no) as (we assume eo = 0) 



-p-^ In 



1 



1 



(9) 



Various thermodynamic quantities of the gas such as the 
chemical potential /i, the internal energy U = J iidN and 
the specific heat C = dU /dT can be obtained analytically 
if we have an analytic expression for (uq). Equation 
shows that fi is nonzero for finite N, in contrast to the 
thermodynamic limit iV — > oo where fi vanishes. 

By using Eq. 0, (no) can be singled out from Eq. (jHl 

as 



which may be rewritten as 
1 



1 



+ 1 e/S^A. 



N - (no) = 



1 



fe>0 



(«o) + l 



(10) 



(11) 



This gives a self-consistent equation for the mean par- 
ticle number (no)- In the limit (no) S> 1, the term 
(no)/((no) + 1) inside the summation may be approxi- 
mated by 1. As in the previous section, we denote the 
summation over k as 



n 



5Z g/3efc 
fc>0 



1 



1' 



(12) 



yielding a quadratic equation for the mean number of 
particles in the ground state 



N ~ (no) = 



n 



1 



whose solution is 

, , 1 
2 



N -n-i + ^/{N -n-i) 



-47V 



(13) 



which exactly coincides with Q, linking the two ap- 
proaches. It is natural to ask about the physical sig- 
nificance of this coincidence. The master equation re- 
sult is derived from a particle-number conserving formal- 
ism, which already has the canonical ensemble property 
built in. However, the master equation is solved approxi- 
mately in Eq. under the condition of large condensate 
particle number. On the other hand, the result of the 
grand canonical analysis, Eq. H13() . is derived by trading 
the chemical potential for the average particle number, 
which corresponds to imposing the total particle number 
constraint on average. The resulting equation is further 
solved self-consistently, also under the large condensate 
number assumption. The coincidence of results in this 
approximation, and the excellent agreement with the nu- 
merical canonical ensemble simulation, indicates that a 



N=200 



canonical 
exact 




FIG. 2: (Color online) The average condensate particle num- 
ber is plotted versus temperature for A'' = 200 particles in 
an isotropic harmonic trap. The approximate analytic result 
given by Eqs. O (or 1131 and 11411 (solid line) is compared 
with the "exact canonical dots" computed numerically with 
good agreement. 



strict accounting of the canonical ensemble constraint is 
unnecessary for the calculation of (no). 

One should mention that the trick used in Eq. Hll|) is 
nontrivial. A straightforward expansion in l/(no) yields 
very poor accuracy. We discuss this in Appendix A. 



IV. ANALYTIC EXPRESSIONS FOR MEAN 
CONDENSATE PARTICLE NUMBER AND 
VARIANCE WITHIN THE GENERALIZED 
GRAND CANONICAL APPROACH 



It is of interest to have an accurate analytic expression 
for the mean number of particles in a mesoscopic conden- 
sate that is valid for all temperatures. Let us specialize 
to the case of an isotropic harmonic trap. Single-particle 
energy levels are then Ck = h,^l{l + m + n), where fl is the 
trap frequency and k = {l,m,n} denotes the quantum 
numbers. In this case, the coefficient Ti, can be evaluated 
approximately by the method shown in the Appendix B. 
The idea of this method is to first convert the triple sum 
into a single sum, and then to approximate the single 
sum as an integral. We obtain 



C(3) 



f- 

\2a 



1 



(In 2 - In a) 



(14) 



where a = phn. The first term has been previously ob- 
tained in Ref. [l^ . and here we find subdominant cor- 
rections in a. Thus, we have an analytical expression for 
(no) given by Eq. (|13|l . where the term N ~ H — I can 
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be expressed as 



N-H-l^Nil- \ — 



where Zk 
that 




f N 



2/3 / rp\ 2 



( ^ 



VC(3) 



1/3 



T 



In 



4 VC(3)y VT, 



T 



In 2 



and the critical temperature Tc = (iV/(^(3))^/'^?iil/fcB has 
been introduced. For large N , the second line of 115|) may 
be neglected, and if we also take (iV — Ti — 1)^ ^ 47V, 
then Eq. ||T3J| reduces to Eq. (8) of Ref. \l2\: 



(no) -Nil 



( ^ 



VC(3) 



2/3 y N 2 

27; 



(16) 



1 - {TIT,) 



in the 



which further reduces to 
thermodynamic limit. 

We now consider fluctuations of the number of parti- 
cles in the condensate, starting with the variance. The 
variance of the condensate particles can be expressed in 
terms of the non-condensate variance using the identity 
no = N — J2k>o''^k- In the grand canonical approach, 



(n|) can be evaluated from 



N=200 



'd^Zk, 



(17) 



<1 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 
T/T 



(1 



It follows from Eq. |[T7|) 



2{nkf 



{Uk) 



(18) 



We then obtain 



(15) 



^2 = ((no - (no))^) w ^[(n^) - (rife) 
fe>o 

~ Y.^{nk)^ + {nk)] 

1 1 



fc>o 



E 

A:>0 



(Ae^*^" - 1)2 AeP'>' - 1 



(19) 



where A = 1 + l/(no)- We have assumed that fluctua- 
tions of Uk (fc 7^ 0) are uncorrelated {{nkUm) = {'nk){nm), 
k ^ m). This is the case provided (no) is much larger 
then its variance, so that particle exchange with the con- 
densate reservoir is the main channel of fluctuations of 
Uk- Near and above Tc, the correlations become sub- 
stantial, which yields the failure of Eq. I|19|l . Taking 
the chemical potential ^ as given in @ and using the 
method described in Appendix B, we find the following 
analytic expression for the variance in the generalized 
grand canonical analysis: 



M2 



In 



A 



7Ae"/2 + 8 
8a(Ae''/2 - 1) ^ 2a2 V^e'^/^ - 1 



(20) 



Ae"/2 - 1 



lnA + dilog(Ae"/2) 



where a = (3hVl and dilog(a;) = dt\n{t)/{l — t). In 
the hmit ksT > hft and (uq) > 1, Eq. yields ^2 ~ 
7r2/6a3 = TT'^T^N/6Tj?C{3) which agrees with Eq. (11) of 
Ref. 

The analytic results are compared with numerical sim- 
ulation in Figs. I2I8I Fig. [21shows an excellent agreement 
between (no) computed analytically from Eqs. (|13|l and 
(|14|l and the exact numerical simulation obtained in the 
canonical ensemble _1.9|. In Fig. Owe plot the variance 
Ano as a function of temperature obtained from Eq. I|2U|) 
(solid line), as well as the "exact canonical dots" obtained 
from numerical computation [l9j. While the analytic re- 
sult is good at low temperature, at temperatures com- 
parable to Tc and above, there is substantial deviation 
from the canonical ensemble result because correlations 
between excited levels are neglected. 



FIG. 3: (Color online) The variance of the condensate distri- 
bution is plotted versus temperature for — 200. The ap- 
proximate analytic result within the generalized grand canon- 
ical approach (solid line) is compared with the "exact dots" 
computed numerically in the canonical ensemble. The an- 
alytic result is in good agreement for low temperature, but 
fails at temperatures comparable to the critical temperature 
and above. The reason for the high-temperature discrepancy 
is that the modified grand canonical theory neglects correla- 
tions between excited levels. 



V. ANALYTIC EXPRESSION FOR ALL 
CUMULANTS VIA THE STOCHASTIC PATH 
INTEGRAL FORMALISM 

While the master equation is a powerful approach, it 
is often the case that a simple analytic solution cannot 
be found. It is therefore of great interest to pursue al- 
ternative treatments that give an approximate analytic 
solution of the master equation, that is asymptotically 
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valid in the physically relevant limit of many particles in 
the condensate, (no) ^ 1. Just such an approach was 
developed in Refs. fill ITfil |. by solving the fluctuation 
statistics problem with a stochastic path integral. The 
stochastic path integral formalism is complimentary to 
the master equation approach as will be shown below. 
However, we stress that it can also be applied in cases 
where it is impossible to even write down a differential 
master equation, as demonstrated in Ref. |16| . 

The calculational details are given in Appendix C, but 
the basic idea is to translate the master equation into the 
stochastic path integral, whose action functional contains 
all rate information, and also imposes local particle con- 
servation. The fluctuation statistics can then be calcu- 
lated in saddle-point approximation by finding the "zero 
energy lines" of the dynamics - the statistical trajectory 
in phase space that is most likely, similar to the instanton 
trajectory of Ref. |2Q|. From this trajectory, the gener- 
ating function may be found as an area in phase space. 

Rather than solving the original master equation 
we skip to the generalization of CNB2 where 

K„„ = (7V-no)(l + 7]), Hr,„=n+{N-no)v, (21) 



which applies also to higher temperatures, and r] is de- 
fined as 

(22) 

fc>o ^ 

We introduce the notation k„ for the cumulants of the 
distribution |2ll |. and define the nonstandard cumulant 
generating function Q{X) as 

i^s=dr'Q{X)U=o- (23) 

According to the calculations in Appendix C, this func- 
tion is given by the solution of the equation 

(Q + l)KQ{e^ - 1) + QHQie-^ - 1) = 0, (24) 

for any K, H . For the special case of (|21() . the result is 



Q(A) = 



-n+{l + rj){N - l)e^ ~Nr]+ ^/4e^{l + r]){-r] + e^(l + ri))N + [U - e^(l + ri){N - I) + iiNf 



-2?7 + 2e-^(l + ?7) 
Applying Eq. the average value (ki = {no)) is given by 



(no) = {l/2)[{N-n-l-ri) + ^(A^ - ^ - r/ - 1)^ + 4iV(l + 77)] , 
which coincides with in the rj = Q limit. The second central moment (k2 = ^2) is given by 

(1 + r;)(r; + n)yJ{N - H - rj - 1)^ + 4N{1 + ??) - (1 + 'n)[T]^ + H{1+H- N)+ f]{l + 2W + A^)] 



M2 



in terms of elementary functions. 



2y'{N -H~f]~l)^ + 47V(1 -I- ry) 
I 



(25) 



(26) 



(27) 



All higher cumulants and central moments may be easily 
computed from l|25|l . Figures comparing these approx- 
imate results with the exact master equation solution 
(shown in Fig. ^) are not shown, simply because for 
A^ — 200, they are indistinguishable. 



VI. HYBRID THEORY - COMBINING CNB3 
WITH THE MASTER EQUATION ANALYSIS 

We now demonstrate how to combine ideas from the 
canonical ensemble quasiparticle formalism of CNB3 0| 
(which works well at low temperature when ^/Jl2 ^ ^o) 
with the physics of the master equation approach, in or- 
der to obtain essentially perfect quantitative agreement 
with the exact numerical solution of the canonical par- 



tition function at all temperatures for the fluctuation 
statistics of the Bose gas. 

Defining the function Fng as the ratio between the 
probabilities to find uq + I and no particles in the ground 

state, 



p — 

^ no 



Pno + 1 
Pno 



(28) 



we note that the canonical ensemble constraint is im- 
posed by F/v = because if all the particles are in the 
condensate, it is impossible to cool further, or, in other 
words, the probability to find A^ -I- 1 particles in the con- 
densate is equal to zero. It is then useful to consider an 
expansion of this function in N — nQ. Rather than Taylor 
expand, a better approach is to approximate this func- 
tion by ratio of two power series and then determine both 
the numerator and denominator coefficients, a procedure 
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FIG. 4: (Color online) The present hybrid theory, we call CNB5, (solid line) is compared with the exact numerical results 
obtained in the canonical ensemble (dots), (left) The average condensate particle number (no) and (right) the variance Ano 
are plotted for A*' = 200 versus temperature. The predictions of CNB3 ^3 (dashed line) are also shown. 



known as a Pade approximation. Pade approximations 
are usually superior to Taylor expansions when functions 
contain poles, because the use of rational functions allows 
them to be well-represented. We approximate 



TP — 

^ nn — 



K 



no — rr 



"0 + 1 



(29) 



where the functions H, K are both polynomials in N—uq, 

Hno = n + r]{N - no) + a{N - nof , 

Kno = il + v){N -no) + a{N -mf, (30) 

and truncate the expansion at second order. Knowledge 
of the function Fng allows the construction of the entire 
distribution, 

Pno = V U ^rn\ = E 11 ^rn\ (31) 



m— no 



no— rn— no 



or 



_ ^ (iV - ?io - 1 + a:i)!(iV - no - 1 + 0:2)! , . 
^""""^ {N~no)l(N-no + {l + v)/c.)l ' ^ ^ 



where 2:1,2 — {rjzt \/ ff' — ^a'H)/2a and C is the normal- 

N 

ization constant determined by ^ = T The func- 

"0=0 

tions _ff , K take the same form as in the master equation, 
but now the coefficients 7i, 77, a are treated as free param- 
eters to be fixed by comparison with CNB3 jl^ at low 
temperatures. 

The further analytic input for the theory is the first 
three moments of the distribution described by the heat- 
ing and cooling coefficients H30|l in the low temperature 
limit. These moments are used to fix the free param- 
eters 7i, 77, a. The calculation is done in Appendix C 
for the complete generating function using the stochastic 
path integral formalism, but here we only reproduce the 
needed first three: 

(no) = N^n, (33) 
/i2 = n{l + ri + an), (34) 
/i3 = -H{l + r] + an){l + 2'q + 'ian). (35) 

Comparison with the CNB3 17] at low temperature al- 
lows us to obtain the parameters 7i , 77, a, 



n 



^=2 



1 



nk 




E 



[nl + nk) 



(36) 



8 



N=200 
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FIG. 5; (Color online) Same as in Fig, 
as well as cumulants k„ 1211. 



01 for the third, fourth, fifth and sixth centered moments, denoted by ((no — i^o))"), 



where fik = {e^'^'' — 1)^^. Knowing these parameters al- 
lows the complete specification of the entire distribution 
in this approximation. To sum up: The theory uses {%) 
the master equation/stochastic path integral formalism 
to determine the first three moments at low tempera- 
ture, and (m) the results of CNB3 (that works well in the 
low temperature limit) to fix the three undefined param- 
eters. The distribution function H32|l together with Eqs. 
(jlffill then give s predictions for all central moments and 
cumulants '2l| at all temperatures. 

The theory is put to the test in Figs. 01 and for the 
first few central moments and cumulants of a thermal 
Bose gas in a harmonic trap, for the mesoscopic case 
N — 200. The hybrid treatment yields perfect agreement 
with the "exact" numerical dots obtained in the canonical 
ensemble. 



VII. CONCLUSIONS 

We have discussed the fiuctuation statistics of an ideal, 
mesoscopic, Bose-Einstein condensate from several dif- 
ferent perspectives. First, we have reviewed the master 
equation approach, and derive an approximate analytic 
solution for the mean condensate particle number. By 
generalizing the grand canonical analysis, the same result 
from the approximate solution of a self-consistent equa- 
tion has been recovered for the mean condensate particle 
number. Improved analytic results are obtained for the 
mean condensate particle number in the case of a 3D har- 
monic trap, that are quite accurate when compared with 
numerical calculation of the canonical partition function 
result. Analogous treatment of the variance of the dis- 
tribution in the generalized grand canonical picture have 
given results that are very accurate below the critical 
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temperature, but substantially deviate at or above the 
critical temperature. The reason for this discrepancy 
is because correlations of excited energy levels were ne- 
glected in the calculation. 

Next, we have presented an (approximate) analytic 
solution for the generating function of the fluctuation 
statistics from the master equation perspective. This is 
done by employing the stochastic path integral formal- 
ism, with the saddle-point approximation, giving results 
that are asymptotically valid in the physically relevant 
limit of many particles in the condensate. The general 
solution is discussed for arbitrary heating and cooling 
coefficients, and specific results are given in terms of ele- 
mentary functions for the Bose- Einstein condensate heat- 
ing and cooling coefficients of CNB2. These results are in 
excellent agreement with exact master equation solution 
when compared numerically. 

A hybrid theory has been put forth that combines 
the master equation/stochastic path integral with the 
approach of CNB3. The theory applies the results of 
CNB3 in the low temperature limit, together with the 
predictions of the master equation/stochastic path inte- 
gral. This is accomplished by preserving the physical 
structure of the master equation, while using the first 
three moments of CNB3 to fix the numerical values of 
the heating and cooling parameters. These predictions 
are then examined for several higher moments (or cumu- 
lants), at all temperatures. The predictions of this theory 
are essentially in perfect agreement with numerical sim- 
ulation of the exact canonical partition function. 

Finally, we briefly discuss how our methods and results 
extend in the presence of weak interactions. As will be 
shown in Ref . [23 , it is easy to generalize the hybrid ap- 
proach to the interacting case and take the first three cen- 
tral moments (in the low-temperature limit) from CNB3. 
From the microscopic master equation point of view, in- 
teractions generally lead to "off-diagonal" transitions in 
the density matrix, demanding a fully coherent treatment 
of the problem. However, it will be demonstrated else- 
where that good agreement with CNB3 (in the applicable 
low temperature limit) may be obtained with only diago- 
nal transitions, where the heating and cooling coefficients 
are determined from Bogolioubov theory. 



APPENDIX A: GRAND CANONICAL 
FORMALISM: SIMPLE EXPANSION WITH 
TRIPLE SUMMATION 

Before calculating the analytic expressions for the av- 
erage and variance of the condensate fluctuations, it is 
instructive to first take another path, different than the 
one presented in Sec. III. Rather than make the step 
we make the simple expansion 



=1 fc>0 ^ 



y\ (no) 



A{T)^—B{T), 

{no) 



where 



k>0 k>0 



(e/3^fe - 1)2- 



(Al) 



(A2) 



The above expansion is in the small parameter 
g/3eiY[(g/3efc _ l)(no)]. Equation I^KT^ has the solution 



(no) = i (iV - A + ^{N-AY+AB) 



(A3) 



For an isotropic harmonic trap, = h^{l -\- m-\- n)^ and 
since e~''('+'"+") < 1 (here a = (iMl) we make the series 
expansion 

°° g-a(/+m+n) 

'^^^^ ~ ^ ^ 2 — g-a('+m+«) 
/,m,n— 1 

00 00 

^ ^ ^ g-(s+l)a(/+m+n) 
/,m,n— 1 s— 



E 
E 



^e-(«+i) 



.ra=l 
1 

a 



E 



E' 



n 3 







rC(3), 



(A4) 



where C,(n) = ^ 4r is the Riemann-zeta function and 

the conversion to integration is based on the assumption 
a < 1. Similarly we find 
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N=200 



canonical 
exact 




is (no) + 1 ^ {e^'^'' — 1)~^, allowing for a better approx- 
imation. 

After this point, there are two more key steps: (1) con- 
vert the triple summation into a single summation, (2) 
make an integral approximation to the single summation. 
We first recall 



n 



1 



^ _ 1)' 



(Bl) 



where for a harmonic trap /Je^ = a{l + rn + n) and we 
have introduced a = Mlf}. The triple summation can be 
reduced to a single summation over s, where s — Z+m+n, 
and weighting this sum with the number of ways W to 
put s quanta into three boxes, W = [s + 2)!/(s!2!) = 
(s + 2)(s + l)/2, 



FIG. 6: (Color online) The average condensate particle num- 
ber is plotted versus temperature for A'' = 200. The ap- 
proximate analytic result 1A6II within the grand canonical 
approach (solid line) is compared with the "exact canonical 
dots". The triple sum, together with the simple expansion 
gives a poor fit. 



Finally, by noting that A{T) = N{T/Tcf, B = 
7V(r/r,)3C(2)/C(3), where T, = [N/Cm^'^'^^/kB is 
critical temperature in the thermodynamic limit, we ob- 
tain 




N 



C(2)£ 
C(3)7V 



.(A6) 



For large N , Eq. (|A6p reproduces the thermodynamic 
limit. However, the triple sum formula is inaccurate for 
small TV especially near Tc- Result ljA6p is plotted in 
Fig. El The fit is poor because N = 200 is not sufhciently 
large. We note that this expansion may be systematically 
improved by keeping more terms in the l/(no) expansion, 
and evaluating the coefficients in similar manner. In the 
next appendix, this analysis is improved. 



APPENDIX B: GRAND CANONICAL 
FORMALISM: IMPROVED EXPANSION WITH 
SINGLE SUMMATION 

The goal of this appendix is to obtain an improved ana- 
lytic expression for both the average and variance of the 
condensate distribution for an isotropic harmonic trap. 
Equation (|Alll is obtained from Eq. IjlOl) by neglecting 
the terms with l/(no)^ and higher under the condition 



(no) > e'='"'/{e' 



1). For the enhanced expansion 



(|ll|l . it is not difficult to see that the validity condition 



n = 



s=l 



(g + 2)(.s + l) 
2(e"''- - 1) 



(B2) 



In order to find an analytical expression for Ti, we in- 
terpret the summation as a Riemann summation, and 
convert it (approximately) to an integral using the mid- 
point rule, /s ~ Io° dsf{s + 1/2). The midpoint 
rule gives the better approximation because it compro- 
mises between the lower and upper summation. Re- 
parameterizing the integral yields 



n 



a/2 

C(3) 



3a; „ 
^ + — + 2 
a"^ a 



a-\-A 



1 

2^ 



'^Aa? 2a2 

2 a 
C(3) 
a3 \2a 



dilog (e2°^ 



1 dx 



_3_ 

16 



,a/2 



- 1 



V + -(ln2-lna) 
/ a 



(B3) 
(B4) 



where dilog(a;) = /^^ ln(t)/(l — t). This derivation gives 
Eq. 1)14(1 . one of our main results. 

Turning to the variance, we proceed in the same man- 
ner to obtain an analytical expression for /i2. The con- 
version from the triple sum to the single sum yields 



= ^ 



(s + 2) (s + 1) 



(^gas _ 2)2 Ae"-^ - 1 



(B5) 

where A = 1 + l/{no). Next, we convert the sum into an 
integral as before to yield 



,/2 (Ae--1)2 2 \a 



l/xVgX^^')^. (B6) 



Then we integrate by parts to find 



15 



1 



Ai2 



8a(Ae°/2 - 1) 2a2 J^,^ Ae^ - 1 



^+3 



dx. (B7) 
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The integral in Eq. (jB7|) can be calculated analytically 
to find 



In 



8a{Ae<^/'^ - 1) 2a^ yAe"/^ - 1 



(B8) 



In A + dilog( Ae 



This derivation gives Eq. (|20|l another main result. 



APPENDIX C: STOCHASTIC PATH INTEGRAL 
SOLUTION OF THE MASTER EQUATION 

The purpose of this appendix is to provide approxi- 
mate expressions for all cumulants of the stationary con- 
densate fluctuations using elementary functions, starting 
from the master equation approach. In order to accom- 
plish this, we employ the stochastic path integral formal- 
ism IHIii. 

Consider a general differential master equation of the 
form 

Pnit) - J2^Wnr,rPrn{t) - WmnPn{t)], (CI) 



where Wnm is a transition rate from m to n, and P„ is the 
probability of occupying state n. Applying this equation 
to our BEC problem, the microscopic rates to different 
states are taken from CNB2 m: 



k(1 - 

K\n 



' 1 - 
no 



no)no, 



1), (C2) 



where ng is the number of particles in the condensate, k 
is a rate constant, H. is given in (|12|) and 77 is given in (|22|l . 
We now express this master equation as a stochastic path 
integral, by going to a continuous representation where 
the discrete number of particles in the ground state no is 
replaced by an effectively continuous variable Qj 

UiQf,Q^,t) = JvQV\exp{S{Q,\)}, (C3) 

S{Q,\) = / dt'[^XQ + H{Q,X)]. (C4) 

The object U is the evolution operator going from one 
particle configuration Qi to another particle configura- 
tion Q/ in time t. It is expressed as a path integral over Q 
and A. The auxiliary variable A is a canonically conjugate 
variable and imposes local particle number conservation. 
In the continuous limit, (and suppressing the overall rate 
constant k) the CNB2 >Q| rates may be expressed as 

WiQ',Q) = il + 7^){N + l-Q')Q'S{Q' -Q-l) (C5) 
+ [H + v{N -Q' - 1)](0' + 1)S{Q' -Q + 1). 



According to the prescription of Ref. the Hamilto- 
nian H{Q, A) of the stochastic path integral is found from 
the equation: 



HiQ,X) = / dQ' 



.{Q'-Q)>^ 



-1 W{Q',Q). (C6) 



For the rates ljC5|l . we find 



H(Q,A) = (l + 77)(iV-0)(Q + l)(e^-l) 

+ [n + ii{N-Q)]Q{e-^-l), (C7) 

or for the general master equation ^ with arbitrary co- 
efficients H,K, we find 



HiQ, A) = Kq{Q + l)(e^ -1) + Hq Qi. 



-A 



1). (C8) 



This result has a simple physical interpretation: On a 
short time scale, the elementary transitions into and out 
of the condensate are Poissonian, witnessed by the gen- 
erators G of Poissonian statistics, G ~ r[exp(±A) — 1] 
(counting an incoming (-I-) or outgoing (— ) boson). 
These boson transitions are described with a rate into 
the condensate Fin = Kq{Q + I), and a rate out of the 
condensate Fout = Hq Q- 

The stochastic path integral (|C3() may be evaluated in 
saddle point approximation, where the large parameter 
of the expansion is (no) 3> 1, the number of particles in 
the condensate. Applying this approximation gives the 
analog of Hamilton's equations of motion. 



Q = dxH, 



A =. -doH. 



(C9) 



To solve the problem of instantaneous particle num- 
ber statistics in this approximation, we generalize the 
method of Ref. 16], following the method of Ref. '20j, by 
first finding the "zero energy lines" , implicitly defined by 
the equation H[Q, A) = 0. For time scales longer than 
the relaxation time k^^, any Q"distributed initial state 
will be projected onto the zero energy lines. The trivial 
zero energy line is given by Aq = and must exist for 
the probability distribution to be normalized. The in- 
stantaneous fluctuation statistics (to leading order) can 
be found by calculating the statistical action (|C4|) along 
the non-trivial zero energy line. This action, S{x), is also 
the generating function of the cumulants of the fluctua- 
tion statistics. On the zero energy line, the Hamiltonian 
vanishes, leaving only the dynamical part of the action. 



Six) 



dt'Xit')Qit') 



QiX)dX, (CIO) 



and we have changed variables from time to phase space 
coordinates. In the case of Bose-Einstein condensation, 
(|C7|) . the nontrivial zero energy line (5(A) is given by 



12 



, -7^ + (1 + ry)(iV - l)e^ -Nf]+ V4e^(l + + e^(l + + {H - e^(l + r;)(iV - 1) + ?7iV)2 

-2r^ + 2e^(l+ry) ' ^^^^^ 

I 



In order to have a generating function, it is unnecessary 
to perform the integral (jClOl) because dS/dx = Q{x)- 
Therefore, all cumulants of the distribution can be found 
from 



Q(A)| 



A=0- 



(C12) 



This result generalizes the discussion of Ref. jl^J for any 
two Poissonian processes in series (equilibrium or not), 
and is easily generalized to arbitrary elementary pro- 
cesses. Equations (|Cllll and (jC12|) recover Eqs. (|23|l and 
(I25|l and are main results. 



In section VI, this same method is used to calculate 
the cumulants where the heating and cooling coefficients 
are given by 



(l + r;)(Af-no) + a(iV-no)' 



Hno = n + viN -na)+aiN -noY, (C13) 



As before, we define Q(A) as the solution of the approxi- 
mate zero energy equation, = Hq/Kq, where the low 
temperature limit allows the approximation Q + 1 ^ Q. 
The solution is 



Q{X) = N- 



?7 - e^(l + 7]) + ^4aH(e^ - 1) -f [77 - (1 + 77)6 
2a{e^ - 1) 

I 



A12 



(CM) 



The first three cumulants are given in Eqs. (|33(l - (|35|l . 
which coincide with the first three central moments. 

We also briefly note that time-averaged fluctuation 
statistics may be easily calculated within the stochas- 
tic path integral formalism. Consider a detector that has 
finite time resolution. The physical quantity that is of in- 
terest is then the condensate particle number, averaged 
over some time window t, 



Qr = (1/r) / dt'no(t') 



(C15) 



where we take r longer than any dynamical time scale, 
for simplicity. Following the method of Ref. , we find 



the distribution P{Qt) is approximately given by the ex- 
pression 



logP(Qr) 



2 

out J ; 



(C16) 



where Tin = Kq^ {Qr + 1), Tout = Hq^ Qr- Interestingly, 
this result is of the same form as the time-averaged elec- 
tron fluctuations on a mesoscopic cavity, out of equilib- 
rium . This similarity of statistics for radically differ- 
ent physical systems originates from the fact that both 
systems can be described as two Poissonian processes in 
series. 
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